{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior Predictive Checks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial you will learn how to run posterior predictive checks in HDDM.\n",
    "\n",
    "A posterior predictive check is a very useful tool when you want to evaluate if your model can reproduce key patterns in your data. Specifically, you can define a summary statistic that describes the pattern you are interested in (e.g. accuracy in your task) and then simulate new data from the posterior of your fitted model. You can the apply the the summary statistic to each of the data sets you simulated from the posterior and see if the model does a good job of reproducing this pattern by comparing the summary statistics from the simulations to the summary statistic caluclated over the model.\n",
    "\n",
    "What is critical is that you do not only get a single summary statistic from the simulations but a whole distribution which captures the uncertainty in our model estimate.\n",
    "\n",
    "Lets do a simple analysis using simulated data. First, import HDDM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import hddm\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate data from known parameters and two conditions (easy and hard)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data, params = hddm.generate.gen_rand_data(\n",
    "    params={\"easy\": {\"v\": 1, \"a\": 2, \"t\": 0.3}, \"hard\": {\"v\": 1, \"a\": 2, \"t\": 0.3}}\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, lets estimate the same model that was used to generate the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No model attribute --> setting up standard HDDM\n",
      "Set model to ddm\n",
      " [-----------------100%-----------------] 1000 of 1000 complete in 12.5 sec"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<pymc.MCMC.MCMC at 0x104066090>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = hddm.HDDM(data, depends_on={\"v\": \"condition\"})\n",
    "m.sample(1000, burn=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll want to simulate data from the model. By default, `post_pred_gen()` will use 500 parameter values from the posterior (i.e. posterior samples) and simulate a different data set for each parameter value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           knode_name stochastic observed   subj        node      tag  \\\n",
      "a                   a       True    False  False           a       ()   \n",
      "v(easy)             v       True    False  False     v(easy)  (easy,)   \n",
      "v(hard)             v       True    False  False     v(hard)  (hard,)   \n",
      "t                   t       True    False  False           t       ()   \n",
      "wfpt(easy)       wfpt      False     True  False  wfpt(easy)  (easy,)   \n",
      "wfpt(hard)       wfpt      False     True  False  wfpt(hard)  (hard,)   \n",
      "\n",
      "                depends hidden   rt response subj_idx condition      mean  \\\n",
      "a                    []  False  NaN      NaN      NaN       NaN   1.90643   \n",
      "v(easy)     [condition]  False  NaN      NaN      NaN      easy   1.06179   \n",
      "v(hard)     [condition]  False  NaN      NaN      NaN      hard  0.943137   \n",
      "t                    []  False  NaN      NaN      NaN       NaN  0.341715   \n",
      "wfpt(easy)  [condition]  False  NaN      NaN      NaN      easy       NaN   \n",
      "wfpt(hard)  [condition]  False  NaN      NaN      NaN      hard       NaN   \n",
      "\n",
      "                  std      2.5q       25q       50q       75q     97.5q  \\\n",
      "a            0.123325   1.70192   1.81497   1.89453   1.97966   2.18115   \n",
      "v(easy)      0.213885  0.672581  0.925245   1.05542    1.1914   1.49813   \n",
      "v(hard)       0.17852  0.605106  0.817572  0.944719   1.06105   1.27941   \n",
      "t           0.0261543  0.278789  0.325973  0.344599  0.362117  0.383125   \n",
      "wfpt(easy)        NaN       NaN       NaN       NaN       NaN       NaN   \n",
      "wfpt(hard)        NaN       NaN       NaN       NaN       NaN       NaN   \n",
      "\n",
      "                mc err  \n",
      "a           0.00666755  \n",
      "v(easy)     0.00718197  \n",
      "v(hard)     0.00616936  \n",
      "t           0.00127485  \n",
      "wfpt(easy)         NaN  \n",
      "wfpt(hard)         NaN  \n"
     ]
    }
   ],
   "source": [
    "print(m.nodes_db)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd5xU1fn/32dmZ3uvsJ2ySO8IiAUFFcWGFUWMGjUaTWKKMck3P00xJt98Y4oxaog1MSr2glgQUUERAel9Fxa299535vz+ODPb2DKzc2dndjjv12tf98695577zOzdz555znOeR0gp0Wg0Gs3wx+RtAzQajUZjDFrQNRqNxk/Qgq7RaDR+ghZ0jUaj8RO0oGs0Go2fEOCtG8fHx8vMzExv3V6j0WiGJdu3by+XUib0ds5rgp6Zmcm2bdu8dXuNRqMZlgghjvd1TrtcNBqNxk/Qgq7RaDR+ghZ0jUaj8RO85kPXaDSnNlarlcrKStra2rxtik9isViIjY3FbDY7fY0WdI1G4xUqKysJDg4mPj4eIYS3zfEppJTU19dTWVlJQkKvAS29ol0uGo3GK7S1tREeHq7FvBeEEISHh7v87UULukaj8RpazPtmMJ+Ndrlo3MZmkzzzxTE+O1zG2MRw7rvwNEID9aOl0Qw1eoSucQubTfK9l3bw0HsH2HiknGe/yOW+13aj8+xrhitFRUVccsklQ3KvxYsXU1VVZVh/WtA1bvHXjw/z3p4iIoIDuPOcMYQGmnlvdxEbj5R72zSNZlD8+c9/5vbbbx+Se61cuZLHH3/csP60oGsGzXu7i3j0k2xMAh67YSY/u2g83zl7DABv7SjwsnUaTf/cf//93cT0V7/6FY888givv/46S5YsAVRo5X333cecOXOYOnUq//znPwGor69n0aJFzJw5kylTpvD2228D0NDQwNKlS5k2bRqTJ09m9erVrF+/nmXLlnXcZ926dVx55ZUAXHbZZbz00kuGvSft6NQMiv2Ftfzk1V0A/PyiCZwzToVWXT49mb98fJgP9xXT1GolJND5GFrNqUvmz97zSL+5f1ja57nly5dz77338t3vfheAV155hSeeeIKYmBiCgoIAePrpp4mKimLr1q20tLSwYMECLrjgAtLS0njzzTeJjIykvLycefPmcdlll/HBBx+QnJzMe++p91NTU0NkZCR33303ZWVlJCQk8Oyzz3LLLbcAEBMTQ0tLCxUVFcTFxbn9fvUIXeMyFfUt3P7vbTS1WblyRgq3nTWq41xmfBjT0qJpaLXyZY52u2h8lxkzZlBaWkphYSG7du0iJiYGi8XSLe77o48+4t///jfTp09n7ty5VFRUcOTIEaSU/OIXv2Dq1KksXryYgoICSkpKmDJlCh9//DH3338/GzduJCoqCiEEK1eu5IUXXqC6uprNmzdz0UUXddwjMTGRwsJCQ96THqFrXKLNauO7//2GguompqVF8/CVU04Kr5o3OpZdedXszKtm0YQkL1mqGU70N5L2JFdffTWvvfYaxcXFLF++nJCQEJqbmzvOSyn5+9//zoUXXtjtuueee46ysjK2b9+OxWIhMzOT5uZmxo0bx/bt21m7di0///nPueCCC3jggQe45ZZbuPTSSwkODuaaa64hIKBTepubmwkJCTHk/WhB17jEr9/dx5ZjlSRGBLFq5SyCLSe7VKanRgOwM696qM3TaFxi+fLl3H777ZSXl/PZZ58RGRlJbm5ux/kLL7yQJ554gvPOOw+LxcLhw4dJSUmhpqaGxMRELBYLGzZs4PhxldG2sLCQ2NhYbrzxRsLDw3nuuecASE5OJjk5mYceeoh169Z19C+lpLi4GKNqQ2hB1zjNi1tO8MJXJwgMMPHPlbNIigzutd20NCXou/KqkVLqxSMan2XSpEnU1dWRkpLCyJEjARgzZgzZ2dmMHTuW2267jdzcXGbOnImUkoSEBN566y1WrFjBpZdeyuzZs5k+fTrjx48HYM+ePdx3332YTCYsFgtPPPFEx71WrFhBWVkZEydO7Di2fft25s2b123E7g5a0DVOsS23kgff2QvA75dNYUZ6TJ9tR0YFkxARRFldC7kVjYyKDxsqMzUal9mzZ0+31/fccw/PPfccDz30ECaTiYcffpiHH374pOs2b9580rHMzMyT3DMONm3adFI45H/+85+OSVkj0JOimgEprG7izhe+oc0quXXBKK6aldpveyEEU1OiABUNo9EMJ5YtW2aYC8TBrFmz2L17NzfeeGO345MnT2bRokWG3UeP0DX90txm5Tv/2U55fQsLxsbxi4vHO3XdmMRw1h8s5WhZvYct1GiM57bbbjO0v+3bt/d63OgFTHqErukTKSU/fW03ewpqSIsN4bHrZxJgdu6RGW13sxwtb/CkiRqNpgt6hH6K0txm5fEN2Rwtb2BaajQr5qWflFDr8U9zeGdXIWGBZp66aQ4xYYFO9z86IRxAj9A1miFEC/opSHOblev/9RU7TqiwwjW7i3j2i2P879VTOSsrASklqz4/yv99eAiAv1w3ndNGRLh0j9EJ9hF6WYOOdNFohggt6KcgL319gh0nqhkZFcwtCzJ5Z1chewtqWfn01ywan0hlY2uH2P/28klcMGmEy/eICwskMjiA2uZ2yupbSIzoPcRRo9EYh/ahn2I0t1l5/NMcAH592STuOHsMb313AfddeBoWs2D9wVJ2nKgmOtTCP26Yycr5mYO6jxCii9tF+9E1mqFgQEEXQqQJITYIIQ4IIfYJIX7QSxshhHhUCJEthNgthJjpGXM17vLZ4TLK6loYPyKC8yeqZfkBZhN3nzuWT+87l0eumcbjK2by+U/PZenUkW7dyxF/fqKi0W27NZqhoms+9Oeee4577rnHsL4//fTTjr7XrFnDgw8+aFjf4NwIvR34sZRyAjAPuFsIMbFHm4uALPvPHcATaHySj/aVAHDptOST/Nop0SFcNSuVi6eMJDLY4va90mJDAThRqQVdM3wwMh+61Wrt89zSpUt55513aGw07u9jQB+6lLIIKLLv1wkhDgApwP4uzS4H/i1VmZqvhBDRQoiR9ms1PkK71cYnB5WgXzjJ80mzMrSga5zlV1Ee6remz1P3338/GRkZHSs1f/WrXxEREcHrr7/OQw891NGusLCQJUuWkJOTw7Jly/jjH/8IwF133cXWrVtpamri6quv5te//jWgVoveeuutfPTRR9xzzz1ER0dz7733Eh8fz8yZnc4LIQQLFy5kzZo1XHvttYa8XZd86EKITGAGsKXHqRQgr8vrfPuxntffIYTYJoTYVlZW5pqlGrfZlV9DVWMbmXGhjLH7tz1JepwS9ONOCHpl8XFy/n03bbteAVvfoxqNxiiWL1/O6tWrO16/8sorzJ49u1s+dICdO3eyevVq9uzZw+rVq8nLU1L3u9/9jm3btrF7924+++wzdu/e3XFNcHAwmzZt4oorruD222/n3XffZePGjRQXF3ezYfbs2WzcuNGw9+R0lIsQIhx4HbhXStlzPXdvMWknFZWUUq4CVgHMnj1bF50cYrblVgIwf0z8kIQRpttH6HkDCPqRkjq+fvKnrOADOPoC1BfDgu973D6ND9HPSNpTdM2HXlZW1ms+dIBFixYRFaW+QUycOJHjx4+TlpbGK6+8wqpVq2hvb6eoqIj9+/czdepUAK677joADh48yKhRo8jKygLgxhtvZNWqVR19G5kLHZwcoQshLCgx/6+U8o1emuQDaV1epwLGWakxhK25qhjtnMy+E2sZSUJ4EEEBJiobWqlrbuuz3T8/P8pM276O1zJ7XZ9tB4WUsOH38NgcePekOX3NKYwjH/rq1at7zYcOdButm81m2tvbOXbsGH/6059Yv349u3fvZunSpd2uCwvrTEjX3+DJyFzo4FyUiwCeBg5IKf/cR7N3gJvs0S7zgBrtP/ctbDbJ9uNqhD4nM3ZI7mkyiY6J0bzKpl7bSCnZffgYE0ydHjt54mtobzHOkLKD8NkfoPwwbH8Omoyrsq4Z3ixfvpyXX36Z1157jauvvppx48Z1y4feF7W1tYSFhREVFUVJSQnvv/9+r+3Gjx/PsWPHyMlRocI964cePnyYyZMnu/0+HDgzQl8ArATOE0LstP9cLIS4Uwhxp73NWuAokA38CzAuH6TGEI6WN1DV2EZSZBCpMcaNCAaic2K091j0nLJ6MhpUbdKjYdM5aEvDZG2Ggm+MMyJnQ/fXeV8b17dmWNMzH3pYWFhHPvT+mDZtGjNmzGDSpEnceuutLFiwoNd2wcHBrFq1iqVLl3LmmWeSkZHR7fyGDRtYutS4ak3ORLlsoncfedc2ErjbKKM0xrOnQK38nJYaPaTL8AcKXdx0pJwZJvXH05J8OlsO5DLelAfHN0HGfGOMOPqp2kYkQ10hHP8SxvWes1pz6tFfPvSbb76Zm2++uePcmjVrOvYd1Yh60nOEv2TJEg4ePHhSu5KSEpqampgyZcqgbe+JXil6irAnX81jT0nxUHhYH6QPIOg786pJFyqUMixlEntso9WJ0gPGGNDeCrmb1P7Cn6ntia+M6Vvjl3giH3pvnDhxgkceecTQPrWgnyLsLVRRBJNTvSPox/tYLXqwuI50UQpAQtpp5MhkAGTZYWMMKD8EbQ0QOxomXqaOFe1UE6UaryN99PdgdD703pgzZw7Tp0/v8/xgPhst6KcANpvsqBw0OXloBT0jru/QxZZ2K9ml9R2CHpI0hqZINUKXFUfAZnPfgFL7V93EiRASA8HR0N4MDeXu961xC4vFQn19vc+KujeRUlJfX4/F4tqKbZ1t8RTgeGUj9S3tJEUGkRARNPAFBpIaowQ9v6oJq01iNnX673NKGwi11REtGsASCmEJJCUlUXo8msT2aqjNh+h09wwos7tuEieobVQaNFdDzQkIT+j7Oo3HiY2NpbKykrq6Om+b4pNYLBZiY12LSNOCfgpwqFiNzieMjBzye4cEmkmMCKK0roWimqYOgQc4UFRLmrCvGI7JBCHISgznaO5IEkW1CjN0V9AdI/QEe+m86DQo2QM1+ZAyy72+NW5hNptPWsSjcQ/tcjkFOFKiqgaNS3KtSIVRdEyM9vCjHyiq7XC3EJMJQGZ8GDk25Uen/Ij7Ny/rIehR9gLXNfnu963R+Bha0E8BDpcqQR+b6Pn8Lb2REadWzR2r6B6L3nVC1CHoKdHBHJX2tL3uCnpbM1QdA2GGeLX0ukPQq/P6vk6jGaZoQT8FOFKifJTeGqGPSVSCnlPaKehSSg4U1ZIs7JOTUSpzRHJ0CMelPRNk9XH3blyZA9Km/lkEBHW7DzVa0DX+hxZ0P6fdauNouRJSb43QHZkdc7oUjC6rb6GioZUUsz0pU6QalSdHh3BCJgIgq9wU9Aq13Jq4sZ3HtKBr/Bgt6H5OXlUTre02kqOCCQ/yzhy44x9JV0E/UKS+NaQH2gU9Qgl6ZLCF6iC7y6X6uHuhi5UOQR/TeUy7XDR+jBZ0P8fhbhnrJXcLqEnRAJOgoLqJplaV6/xAkYq8SUSlJCCisxB1bHQMZTISYW2FOjdyvFUetXc4uvNYeBIIEzRVqlWkGo0foQXdzzlinxDN8pK7BcBiNpERF4qUcMzu/tlbUIPARmR7hWoU3inoydEh5NndLlTlDv7GFb0IuskEofFqv1EvLtL4F1rQ/ZzOCVHvCTp0ul0O2+3ZcaKaOOowyXYIiQVLcEfb5OjgTkF3Z2K0N5cLQLi97/rSwfet0fggWtD9nCMdIYvec7kATE2NBuCbE1WU1bVQUN1ERg//uYOuE6ODHqG3Nih3jckCkandz4XZF7M06DKIGv9CC7ofY7VJsr0cg+5gdoaqkrQ1t4qdecpvPjfe7sPu4j8HSDFC0CuPqW1MJph7TAY7BF2P0DV+hhZ0P6agqomWdhtJkUFEhbiW5MdopqVFYzELDhbXsuGQEtKpUfYqRj1G6ClG+ND7crdAp8tFj9A1foYWdD/G4a/O8rK7BSDYYmZyShRSwotbTgAwKdKeCqDHCD05OoQ8m30UPdhYdEcMetcJUQfa5aLxU7Sg+zEdES5enhB1sHRK50h8XFI4qY5FRT0EPTEiiFJTPG3SDPXF0Np7LvV+6S1k0YGeFNX4KVrQ/Zgjpb4zQge4dcEoLpiolvV/77wsRH2xOtHD5RJgNpEYGUq+tIcXVp9w/WYOQe/N5RLmcLloQdf4Fzp9rh+T7WMjdJNJ8MSNszhW3qAmab+yLxrqIehgD11sSGQUJcqPnjjetZv153Jx5EHXRS40foYeofspNpvsSJs7NsE3BB3AbBKdETd19hF6ZG+CHjL4WPTWBuWqMQd25m7pio5y0fgpWtD9lILqJprarMSHBxETFuhtc07G2mYXVNHpAulCt1h0Rwiis1Rkq21MJpjMJ593CHpjuTFl7jQaH0ELup/icLd4e4Von9SXAlJNUPaME0eFLh6T9snSChfzopfsU9vEib2fN1sgOEql1m2udq1vjcaH0YLup3ROiPqooDvcLT0iXBwkRweTIx2Viw671rdD0JMm990mNE5tGytc61uj8WG0oPspHf5zL2ZZ7Je6vidEobPQRTtmlerWldDFDkGf1HcbR4IuPTGq8SO0oPsph30gy2K/dAh6XyP0ENoJsFcvkp1+cWdwRtDDdMZFjf+hBd0PkVKSXTJcXC69j9Ajgy1EBAVwxJaiDjjrdqkvU/HlgREQnd53u9BYtdUuF40foQXdDymqaaah1UpcWCBx4UHeNqd3agvVto8ROsDI6GByOgpGOynox79Q2+TpIETf7bTLReOHaEH3Q474SIbFfqktUNueqW27kBwdwiGbPY68cKdz/WZ/rLZjzuu/XcekaKVz/Wo0wwAt6H6Io6iFr6wQ7RXHCD0qpc8mydEhbLOdpl6c+Aps1v77lBKy16v9sYv7b6t96Bo/RAu6H+KIcPGVHC4nIWWXEXpyn81SokMoJJ7qwBHQUgOl+/vvN38b1BWquqEjpvTftqMMnfaha/wHLeh+SEcMuq+O0Juroa0RAsMhKLLPZinRIQAcCrKL8/Ev++/309+r7bTl/fvPodPlon3oGj9CJ+fyM6SUXQpDe3GEXlukkmplzO/lnN3dEpnSr/BmxIUCsNk2gbmsg31vwtzvqJPWNtjzmhq1B0VA2UHIWa/+QSy4d2D7wvTCIo3/oQXdzyita6GuuZ3oUAvx4V7K4bL/HXjzO2oUftXTMOXq7udrBna3AGTEhQHwUv1MfhAShTixWfnSQ2LhzTugcEf3C0wWWPpIZ0hif+iVoho/RAu6n3G4S/y5GMjt4AmkhI9+qcQc4JPfwoTLIKDLP5cO/3nfE6IAMaEWIoIDKGmGpvm3EvrVX+CZC0GYQVpVJsWZ34K2BjAHwYRLBvadOwgMV9e0NapVqIGhg3izGo1vMaCgCyGeAS4BSqWUJyXHEEIsBN4GHCnx3pBS/sZIIzXO07Hk31vulqJdKt1taDyERKsVnofWwqQrOts4BL2fCBcAIQQZcaHsLajl0NhvM6N8jwpLlDaYvgKW/F4l2RoMQqhIl9oCNUrXgq7xA5yZFH0OWDJAm41Syun2Hy3mXuSIt7Ms7n9bbSctgxkr1f7hD7u3cVQgGmCEDp1ul9w6ASteg7s2w0+OwBWPD17MHXSsFtUToxr/YEBBl1J+DujVF8OEbG+XnXMs7JlwCYy70H5sXfe84+X2dLjxWQN2lxGrRs655Y1qVJ00sbPikLvo0EWNn2FU2OJ8IcQuIcT7Qog+MyIJIe4QQmwTQmwrK9MV141GSsnhEi+WnWtrVlEnwgQpsyFhvPJzN5RB0Q6HkZ2JtuLHDdjlqHg1QnfkdzeUjtBFLega/8AIQf8GyJBSTgP+DrzVV0Mp5Sop5Wwp5eyEBINGWZoOyupbqGlqIyI4gMQIL+RwKdkLtnaIPw2CwtWIOusCde7IOrWtL4GWWgiO7hTUfpicotwqewpqjLdXrxbV+BluC7qUslZKWW/fXwtYhBDxblumcZnsks6UuV6JcCn4Rm1TZnYec7hdHH70ru4WJ2zMSgwnKMDEicpGqhtbDTQW90MXGyrg0PtQV2KcTRqNG7gt6EKIEcKuHkKI0+196u+wXqBzQtRL/vNCu6Anz+g8lnkWBASrc/WlneXknHC3AASYTUxKVqtJDR+lu7NatKkanl0CLy2Hv06BnE+MtU2jGQQDCroQ4iVgM3CaECJfCPFtIcSdQog77U2uBvYKIXYBjwLLpZTScyZr+sKx5N9rWRaL96rtyOmdxwJDlagDHFwDZYfUftxYp7udmhoNwO58gwU9zI1J0Q9+1pnS19oCb96pMzdqvI4zUS7XSylHSiktUspUKeXTUsonpZRP2s8/JqWcJKWcJqWcJ6UcIOGGxlN0JOXyxgjdZu0cfSf0GH1PuUZtd7wAB99T+2mnO9319DQl6JuOGOzrHqzLpbkW9r4BCLhnG6TPV3MDX68y1j6NxkV0ci4/4og3y87V5EF7M4SPODk+fMKlKsdKwXbVLioN0s9wuuvzJiQSGGDiq2MVFNc0dzvX0m7FZhvkF8LBhi0eXKNG5RkL1FzAuf+jjm99GtpbBmeLRmMAWtD9hIr6FiobWgkPCmBkVPDQG1Dex+gclNtl5k2dr6deCybnH73IYAuLxiciJbz+TT4Ax8obuO35bUx84ENmPrSOVZ/nuG7zYH3o+95UW0eOmswzIXGSKn13cI3rdmg0BqFzufgJXasUeSXCxeEb72uyc/Gvld+8YBvMv8fl7q+Zncr7e4v587rD7C2oYcOhUprb1GKl6sY2Hl57EIvZxC0LRjnfaWgsIKCpSrmMTOaBr7FZ4fhmte+I4BECZq5UfvV9b8Lkq1x7cxqNQegRup/gVXcLdE4Q9iXo5gCYfQtc/g/nsiH24NzTElkxNx2rTfL+3mKa22wsm5HCtl8u5pFrpgHwpw8PuRbaaDJDSAwglag7Q8k+aK2D6Izu2SInXq62R9ZBS53zNmg0BqJH6H6C18vOlbsWjugqQgh+fdkk5mTGcryikbPGxTMzPQaAq2al8tbOAjYeKeeZL3L50fku2BAaB02Vyu0S5sTyiRNfqW16jzzvkcnq2InNKua+Z8pgjWYI0CN0P8HrZefKB3C5GECA2cQVM1L4weKsDjF38N2FKgzytW15uBQ162ro4gm7uyV97snnTrtYbR2rYjWaIUYLup/Q1Yc+5DRUKEG0hA1YtMJTzB0VS2JEEIU1zewvqnX+wo7QRScnRh2Lp1LnnHzOkeYg++Puycg0miFCC7ofUNXQSnl9C6GB5o46nENKhWvL+T2BySRYNCEJgI/3lzp/oSux6M21qqyeyaLy1fQk4TQVktlY3pmMTKMZQrSg+wHZZZ2jc5PJC4LqmBBN6EXkhpDzJyYC8MmhQQi6MxkXS/erbcL47hWYHAgBWeer/SMfO2+DRmMQWtD9AEfZOa8t+e8IWRw4v7knmZMZixCwv7CG5jarcxe5knGxeI/ajjipcFcnY+2Cnq396JqhRwu6H+D9CVHPRrg4S0SwhbEJ4bRZJQec9aO7slq0ZJ/aJvUj6KPOBnMg5G/TedY1Q44WdD8g22di0L3rcoHOvC+78qqdu8CV1aIl9uRj/Y3Qg8Ih4wxAQs5652zQaAxCC7of4Miy6JW0uW3Nqii0MEOsC6s0PcQ0h6A7m5kxzMlJUZsNSuw+9KQp/bcdu1hts7Wga4YWLejDnJqmNkpqWwi2mEiJ8UKES2UOSBvEZEKAF6ok9WBaR6pdF0foAwl61TFoa4CIkZ3/BPrCIeg5n+jwRc2QogV9mOMoCj0mIRzzKRzh4iArKRwh4HhFI21WJ8S0qw+9vwVJjgnR/vznDhLGQ0SyStZVsmfg9hqNQWhBH+YcKfG2/7xLDLoPEGwxkxoTQrtNcryiceALAkPVgqj25v5zsDjjP3cgBIxdpPa120UzhGhBH+Z0JOXyVtm5gbIseoExCeqfm2OyeEAiR6ptXVHfbRzVmJwZoYMWdI1X0II+zPF+lkXfE/SxdkHPKXNW0O3pCmoL+m7TMUIfYELUweiFIEyQ95VaYarRDAFa0Ic5nVkWvTBCt7ZDmW/50KFzgVWOsyP0CIegF/Z+vrFSVVoKCHG+FmpIDKTMBls7HPvcuWs0GjfRgj6MqWtuo6immcAAE2leiXA5qkqxRaWdXHbOi4yxC3q2yyP0PgTdMTpPmuhcEQwHjmRdh953/hqNxg20oA9jHD7i0fFhBJi98Kt05DZJnDj09+6HUfFhAOSWNzh3wUAul44l/066WxxMuFRtD72nvs1oNB5GC/owxuE/98qCIugU9CTfEvS4sEDCAs3UNrc7V8EoMkVt+xqhD1bQE06DuCxVDen4F65dq9EMAi3owxivL/l35DbxsRG6EIL0ODVKdyp0cSCXiyPCZcRUVw3pHKUfeNe1azWaQaAFfRhz2Ntl5xy+ZR8TdICM2FAAjlc6I+iOEXovLpf2Vig7CIjBvU+HoB9co1eNajyOFvRhjGNR0VhvZFlsrFTFHgJCfCrCxUFGnF3QnfGjh8apDIlNVdDao33ZQbC1QexolXjLVZJnqEnjuiIo2O769RqNC2hBH6Y0tLRTUN2ExSw6xGtIcZRiGzkNzJahv/8ApMe5MEI3mSDGnljMsfLVwWD95w66ul12vzy4PjQaJ9GCPkxxLJoZHR+OxRsRLgV2QU+ZNfT3doKMWOVDP+GMDx0gwb4wypGbxoG7gg4w40a13bW6//QCGo2baEEfphx2uFu85T93uA9SZnrn/gPQ4XKpdDJ00ZHL3ZHKwEHJICdEu5I0CdLPgNY62Pni4PvRaAZAC/owxZED3SsRLjYb5G1R+z46Qh8ZFUyASVBS2+JcObqE8Wpb3kXQrW2d30RGTnPPoPnfVduNj0Crk98aNBoX0YI+TMn2Ztm54t1qAjEqXeVB90ECzCZS7atnTzjjR3e4XMq6uFwKd6oc6HFjISLJPYPGX6ImSOtLYMuT7vWl0fSBFvRhSmeWRS+M0I9+qrajz1GTfj6KS7HocVmAUAU72u2LkY5vUtuMBe4bIwQselDtf/FX9Q9RozEYLejDkKZWK3lVjQSYBJl20RpSOgR94dDf29pEbyIAACAASURBVAU6YtErnPCjB4aqjJG2dsjfqo7l2ld3Zp5pjEFjzlVFpJtr4ItHjelTo+mCFvRhSE5ZPVJCZnwYgQEe+BVWn4Ca/N7PNVZC7iaVGnb0QuPvbSCOiVGnXC7QpRboOjWCPvY5ICDzLOOMcozStzwJdcXG9avRoAV9WOKxCVEpYc2P4K9T4C+T4d17T04qte8NtdBm9LkQFm/s/Q0mvWOE7qygO4pSfAx7XlOZJEef01kAwwhSZyt/elsjbP6Hcf1qNGhBH5Y4criMNVrQj3wE255WqyZNAbD9WXj3B521NqXsDLubep2x9/YAGXZ3lNMj9IwFYAlVsedrf6KOTb/ReMMW3Ku2u17WWRg1hjKgoAshnhFClAoh9vZxXgghHhVCZAshdgshfDMw2Y/wiKBLCRseVvuLfwU3v6eW9e98AT76pTp/6H0Vfx4SA+OXGndvD+EYoedXNWK19VMA2oElGM6+r/N1wniYcInxhqXOVpEzDaWQ84nx/WtOWZwZoT8HLOnn/EVAlv3nDuAJ983S9IdD0B21Mw2h9AAU7VR5TWbfCulzYfkLYLLA5sfg5RtgjX1kec79g8trMsSEBJpJjAiizSoprG5y7qIzfwjn/AymXAu3vA8WDxQOEQKmLVf7+98yvn/NKcuAgi6l/Byo7KfJ5cC/peIrIFoIYaDTUdOVNquN4xWNCGGwoB/5UG3HLekUsbGL4ap/AQIOrVUx1GlzYfa3jbuvh3F5YlQIOPfn6n2HxnrOsKwL1fboZ50uLY3GTQIM6CMFyOvyOt9+7KQS6kKIO1CjeNLT0w249anH8YpG2m2S1JgQQgJdKIc2EIc/UltH2TQHk5ap8nL734bIVDjjexAQaNx9PUx6bBhbc6s4XtHIAifLgQ4JSZMhNB5q86EiB+J9yTjNcMUIQe9tZUmvQw4p5SpgFcDs2bP1sGQQeMR/3lyrlvKbAlSsdE/GnKd+hiEu53QZKkwmFUGz93U4ukELusYQjIhyyQfSurxOBfoo/aJxF0eWRUPdLYXfgLSqfCU+VOzZCDpcLs6GLg4ljvh2R14cjcZNjBD0d4Cb7NEu84AaKeVJ7haNMXhkhJ6/TW19NNGWO7gciz6UOBJ+OVL0ajRuMqDLRQjxErAQiBdC5AMPAhYAKeWTwFrgYiAbaARu8ZSxms4RuqGC3pHbfLZxffoIXWPRpZQIX8o9kzgRhFnlYG9tVOkHNBo3GFDQpZTXD3BeAncbZpGmT6SU5DhG6Ea5XKSEAvsIPdX/BD0m1EJEcAB1ze1UNrQSFx7kbZM6sQSr8n2l+9WPH37+mqFFrxQdRhTVNNPQaiU2LJCYMIMiTWryVThicLSqm+lnCCG6TIz6oNvFUTijeLd37dD4BVrQhxHZRo/OoUvloVk+nQrXHVwuRzeUOErbFfe6EFujcQkt6MOIjggXQ/3n/utucdBRMNoXBT3eXlijItu7dmj8Ai3owwjPRLh0GaH7KS7lRR9q4saobeVR79qh8Qu0oA8jDBd0a7vK3wJ+LejpvuxDj05XkS41+dDW7G1rNMMcLejDiM5FRQZVKSo7oPJyR2f4fG5zd3CELh4rb0D6Wt4Us0WJOhKqcr1tjWaYowV9mFDd2Ep5fSshFjPJUQZlAMz3f/85QHJUMNGhFiobWilwNuviUNLhdsnxrh2aYY8W9GFC54RoGCaTQdEoHREu/i3oQgimpUYDsDOvuuO4lJLnv8zllme/5s/rDntv9B5rF/QKLega99CCPkzweMiinzM9zS7oJzoF/ZVteTz4zj42HCrj0fVH2HKsvyzRHkSP0DUGoQV9mGD4hGhLnSpqYQqAkVON6dOHmZ6uBH1XvhL0wuomfvXO/m5t/rHBS6GDMZlqW33CO/fX+A1a0IcJOWUq5M6wLIv5WwGpFrZ4oiqPjzEjLRqTUC6Xktpmfv/+QZrarFw8ZQQ7HzifEIuZjUfKKa7xQqRJtL02gBZ0jZtoQR8mGD5CP75ZbdPPMKY/Hyc6NJALJ42gzSqZ+/B63t1VSFCAiV9cPIHo0EBOH6WqE32d6wW3S5Q9+3R1HthsQ39/jd+gBX0Y0NxmJa+qEbNJdITguc0Ju6BnnBqCDnDbWd1z1Txw6URSY1SMukPQt3rDjx4Urmq5WltU4WiNZpAYUbFI42GOljUgJWTEhxIYYMD/4PZWu8sFSJ/vfn/DhFkZMdy/ZDzbj1cxMyOaG07vLIPYMUL31sRodDo0Vii3S8QI79igGfZoQR8GdORAN8p/fvwLaG+GhAkQFmdMn8OEuxaO6fX41NQoAgNMHCqpo6apjagQy9AaFp0OhTuUoKedPrT31vgN2uUyDHD4zw1LyrX/bbWdcIkx/fkBQQFmxo+IAOBgUe3QG9AxMXp86O+t8Ru0oA8DjpTWAZBlhKDbrHDgXbU/8XL3+/MjJoyIBOCAVwQ9Q211pIvGDbSgDwOOlKgR+rikCPc72/0KNJar1YlJk93vz4+YMFJ9vgeK6ob+5jp0UWMAWtB9nDarjWPlDQhhQAx6fSms/7XaP/s+vy1oMVjGj1Qj9IPF3nS5aEHXDB4t6D5ObnkD7TZJWkwoIYHmwXeUswGePh/qitRS/6nXGWekn+BwuRwqqcNqG+K8LjoWXWMAWtB9nCP2CdFB+89b6mH1jfCfK1R61pHT4PrVYNK/+p5EhVpIjgqmuU19KxpSdCy6xgD0X7WPc7jEPiE6GP+5zQZv3K4mQS1hsOgBuPUjCE8w2Er/YYJ2u2iGMVrQfRy3Ruj73oBDayE4Gr7zGZz1Y7AEG2yhf+EQdO9EumhB17iHFnQf54h9hO5yhIuU8OXf1f6iByA+y2DL/JPxIx2x6N6MdNGx6JrBoQXdh3FEuIAqbOES+dtUvdDQOJh+gwes80+8O0LXsega99CC7sMcr2igzSpJiw0hNNDFLA2H1qrtlGtPifS4RpEZF0awxURhTTPVja1De3PtctG4iRZ0H8axoCgrcRATokc+UttxFxhokf9jNglOs7u3DhYPsdtFC7rGTbSg+zCHHYKe5OKEaE0BlOwFSyhkLPCAZf6N19wuOhZd4yZa0H2YfYU1QOeCF6fJ3ai2o86GgCCDrfJ/HEm6+hL0VZ/nsOzxL/hgb5GxN9ax6Bo30YLuw+wpUII+JTXKtQsduc7T5hps0alBZyz6yS6XtXuKeHjtQXacqObOF75hc06FsTfXbheNG2hB91FK65opqmkmPCiAUa5WKXIIeuoc4w07BXDkdDlUXEe7tdP10W618dCa7oWlX/zaYOHVgq5xAy3oPspe++h8UnIkJpMLSbRaG6F4LwgTJM/wkHX+TVSIhZToEFrabR0LuwDW7S+hsKaZ0fFhbLr/XISAD/cVGxsNo2PRNW6gBd1H2Z2vBH2qq+6Wop0grZA4SflkNYPCUZJuy9FOl8pzX+YCcNP8DFJjQjljTByt7TY2HDLQ361j0TVuoAXdC0gpeXVbHj9cvZM9duHuyZajqrbljPQY1zrvcLfMdsfEU555o5Wgf2X/PRwoqmXLsUrCAs1cNSsVgLOyVE4cQ+uQapeLxg10TVEv8Lv3DvDUpmMAvL2zgKe+NZvzxid1nG9us7L9RBVCwPzRLtb81P5zQ5g/Oh6ALccqsNkkz9tH51fPSiUiWNUb7RjFa0HX+AhOjdCFEEuEEIeEENlCiJ/1cn6hEKJGCLHT/vOA8ab6BweLa3nmCyXmYxLCsEn48Su7KK1t7mizLbeK1nYbE0dGEhMW6HznUkKeFnQjSIsNISU6hKrGNp79MpfXtucjBNx0RmZHmykpUYRYzBwta6CsrsWYG3cIuo5F17jOgIIuhDAD/wAuAiYC1wshJvbSdKOUcrr95zcG2+k3PPLRYWwSvjU/g3U/PIezsuKpamzjt+8d6Gjz+ZEyAM4Y4+LovLYA6oshOArixhpp9imHEIKb7eL92zX7abdJrp2V1q1qlMVsYmZGNADbjxs0Sg8Mg9B4HYuuGRTOjNBPB7KllEellK3Ay4CuLjwIimuaWX+gBItZcM95WZhMgoeXTSHYYuLdXYV8eqiU5jYrr23PB+DCSSNcu4HD3ZIyWxewMICV8zNIiVZ5cFKiQ/jRBeNOajM9TQn67j7mQgaFdrtoBokzf/UpQF6X1/n2Yz2ZL4TYJYR4XwgxqbeOhBB3CCG2CSG2lZWVDcLc4c1r2/OwSTh/YhIJEWoFZ1psKD9YpITiJ6/u5jdr9lPZ0MrEkZHMynB1QnSb2mp3iyEEW8ys/s48nr/1dNb/+BySIk/OJT8lRQm6YxGYIWhB1wwSZwS9tyDongUXvwEypJTTgL8Db/XWkZRylZRytpRydkLCqVU1x2aTrN6m/i9eNye927k7zh7N3FGxlNe38OIW9Uf87TNHIVwt4qwnRA0nNSaUc8YlEGzpvZ6rI6x0T0ENUhpUh1THomsGiTOCng+kdXmdChR2bSClrJVS1tv31wIWIUS8YVb6AZuPVpBX2URKdAhnju3+0ZhNglUrZ7N8ThqJEUH85vJJXDmzty9B/dDeCoU71X7KTIOs1gzEyKhg4sICqW5sI7+qyZhOnR2hN1XBR/8P/nkO/Os8WPcg1BqcX0YzrHAmbHErkCWEGAUUAMuBbhUThBAjgBIppRRCnI76R2Fwkovhzctb1ej8mtmpmHtZ+RkVauEPV00d/A1K9qiJtLgsCI0dfD8alxBCMCU1ik8PlbE7v4a02FD3O40ZpbblR/puU54Nzy1Vk+AOCrbDtmfg0r/B5Cvdt0Mz7BhwhC6lbAfuAT4EDgCvSCn3CSHuFELcaW92NbBXCLELeBRYLg37/jn8qWpo5cO9xQgB18xOG/iCwZD3tdpqd8uQMzWl0+1iCCMmq23xXhWK2pOmKnhhmRLz1Dlw0zuw8i3IuhBaauG1W2D7c8bYohlWOLWwyO5GWdvj2JNd9h8DHjPWNP/hzR0FtFptnD0uoSNqwnCOf6m26fM807+mT6akOiZGq43pMDxJhS42lkNNXqcLxsH63yh3zMjpcNPbKtQRYPRC+OJv8PGDsOaHaqQ/+hxjbNIMC3Rsm4dpt9p44Ss1uXWdp0bnUsKJzWo/4wzP3EPTJ1McI/R8gyZGheg+Su9K8R7Y9iyYAmDZk51i7rjuzHvhrB+DtMHrt0GTQf9kNMMCLege5q2dhRwtbyA9NpQLJiUNfMFgqMiBhjIIS9ALirxAUmQQCRFB1Da3c6Ky0aBO7YJe0kPQP/8/QMKc2yBxQu/Xnvs/kDZPLUz65LfG2KMZFuhcLgayt6CG17bnU9/SzrTUKAIDTDxkXwF67+IsLGYP/f/M/Vxt0+epUZpmSBFCMDUlivUHS9mdX0OGq/nre2OEfYK84JvOY6UHYf87YA6EBT/o+1qTGS75Czx5ppoknfddiBvjvk0an0cLukHszKtm+arNNLep/BuO1Z4ASyaN4PLpLoYhusKhD9R27Pmeu4emXybbBX1PQQ2XTkt2v8PMM9X26KfQ2qBcKxv/BEiYsRIiB7hH0kSYfgPs+A9sfASueNx9mzQ+jxZ0A2i32rj7v9/Q3Gbj4ikjmDsqjp151VQ3tjI7M5Y7zh7da6iiIbQ2qD96BIxb4pl7aAakY4GRUSkAolJUCoeCbZC9XrnS9r6ufOdn3utcH2f9GHa+CLtehvP+H0SONMY2jc+iBd0A1h8spaC6iVHxYfxt+QwsZhPfGqqbH/5QxZ+nzoEID/noNQPimBjdW1CDzSZdqzLVFxMuVYK+4WEl5NIGs289OeqlL2JHwfiL4cC7sOMFOOc+923S+DR6UtQAHFEsK+amG+8nlxLKDkPhDmhvOfncl4+q/SnXGntfjUskRgYzIjKYupZ2cisajOl06nUQEgtlB9TCschUWPSga33MvlVttz+n0/GeAmhBd5Pqxla+yC7HYhZcba9kYxj52+CJBfCPObBqIfxvJrx/f+fy7t2rldCHJcDMlcbeW+MyU+xul515BoUKRo6ElW9C7GjIugBufB2CI13rY9RCiEqD2nzI/9oYuzQ+ixZ0N/kypwKbhFkZMUSHulCMYiCOb4bnL4XSfWqRSVwWtDXClifhb9PgqcXw9t2q7Tn3g8VDC5Y0TjPXXsFoU3a5cZ0mT4fv74AVr0LieNevN5lgwmVqf/87xtml8Um0oLvJRnsxCkd9SUNorlWLQtoaYdr18KP98L1tcNeXMPFy5TPP3wq2djj7PhWTrPE6jmdg05Fy4zIvGsFEe/mCA+/0nkpA4zfoSVE3+fywGo2dbaSgf/a/6ity8gy47DEw239NSZPg2n9DXYlacJI4YeDwNc2QMS4pnMSIIErrWjhcUs9pIyJ6bSelpLa5ncjgANdTJA+G1DkQMVKlESj8BlJmef6eGq+gBd0NimuaKahuIiIogEnJLvo2+6KxUi0GAbjkr51i3pWIJB3R4oMIITgrK4HXv8nnnV0F3Deiu4skt7yBv358mHX7S2hotRJiMbN4YhL3Ls7qVtrOcEwmGH8JbP2XcrtoQfdbtMvFDXbmVQEwLS3amDA1gK1PK1fL2MXKf6oZVlx/usrX8+KWEzS1WgE1Il+99QQXP7qRt3YW0tBqJdhioqnNyru7Cln66Ebe3VXYX7fuM9HuR9duF79Gj9DdYMcJFc0wIz3auE73vq62c+8yrk/NkDErI4ZpqVHsyq/hN2v2c+c5o3l47QE+3FcCwCVTR3L/kvGkxYaSX9XIIx8d5s0dBfzg5R0EmAQXTfHQ4p/0M9TkeuVRKNnXmfxL41foEbob7LCHpzkKBbtNRY6KOQ6OglFnG9OnZkgRQvCLiycQaDbx0tcnOOf/PuXDfSVEBAXwl+um8ffrZ3QUwUiNCeUv103nB4uysEn48au7OFxS5xnDzAEw4RK1v//t/tvWlajVqW3NnrFF4zG0oA+SdqutY5m3YYJ+4F21HbcEAgwMgdQMKXNHx/H4ipmMjg8jxGLmsmnJrP3BWSybkdrrJOi9i7O4fHoyja1Wfrh6J+1WDy0AckS79CfoRbtVUq8XrlThsRU5rt1DSrC2Dd5GjVtol8sgOVRSR1OblfTYUOLCg4zp1CHo4y8xpj+N11g8MYnFE52buBZC8PCyKWzLrWJfYS1PbzrGd87xQHbEzLMgJAbKD6nMjT3j2ttbYPWNKu0uqIpIa38CN77hXBbPrU/DJw9BU6WaA1ryB4jPMv59aPpEj9AHiWM1oGH+89pClbcjIATGLjKmT82wISwogN8tU37txzZkU93YavxNzBYYv1Tt9zZK3/JPqD4OCRPgRwcgOBpyPoGc9QP3/dWT8N6PlJgDZH8MT5/fPf2vxuNoQR8kjglRw9wtB99T27GLuleh0ZwyLDwtkTPHxlPX3M6Tnx31zE0mXqG2PQW9oQI+/5Pav+Ahtb7hjO+p11uf6b/P6jxVFg9UgeqfHoNxF6napy9eBzX5/V+vMQwt6IOkc4Qeow7sfhX+eQ48djp8+XfXEyFpd4sGuO/C0wB49otjlNR6YFJy1Dlq0r10HxTu7Dz++R+hpQbGnAdZi9WxGStVlsfD70NNQd99fvoHaGtQPvpZN0NoLFz3HzWx31AKr90K1nbj34vmJLSgD4KapjayS+sJNJuYMDIC9r0Jb9wGRTuVf/KjX8I79zgv6o2VkLtJ/fGMu9Czxmt8mmlp0SyZNIKWdhv/2JBt/A0CApVQA3xlL3pRng1bnwJhUqNzBxFJykUjbbDrpd77qytWSeIQ3TNBmi1wzfMQkQx5W9Q/DI3H0YI+CHbnq9H5pJRIglqq4W37V9OFv1APsSUUdv4XtjzhXIeHPwBpVVVqQmM9ZLVmuPCjC8YB8Mq2PKoaPOBLn/sdJd57X4f87fD+fSov0PQVKr1EV6avUNtdL/e+IOnrVWBrU8Lfs8xdaCxc+U9AqFqouV8Y/1403dCCPgi6+c83/Rla62DMIjjnpzDpCrjyX6rhugfV6GcgDqxRW+1u0QDjkiI4Z1wCzW02Xvz6hPE3iE6H2d9WIv7UeWriMzgazvvlyW3HLIKwRKg4AgXbu59rbVCRLQBnfL/3e406W1VOkjZ443b1bVTjMbSgDwKH//z0EaLzgV70QGdo14RLYMaNauTy0f/031lrQ2cUgSMCQXPK8+0zRwHw7825tLZ7IC79godUGCMon/qKVyFixMntzAEw1V48ZeeL3c/tfBGaq1Xyr/S5fd9r4c9Um9oCePf7OvWAB9GC7iJSyg5Bn1/zAbQ3qYmknnlXFj0IgRHKnZLdT9hX9sfQ3qweeJ05UWPnrKx4xiWFU1Lbwnt7PJDnxRIM33oXfrhPhSimnd5322nXq+3e1zurZtmssPkxtT//nv7vZbbAVU9BUKSa/N/+rPv2a3pFC7qL5JQ1UNnQSkKYhah9z6uDc24/uWF4YmcNxw9/0ffqOUfRAe1u0XRBCMGtC9Qo/ZlNuZ7Jry4ERKUOHCY7YjIkTVGj8X1vqmMH34OqXIjOULVPByImEy75i9r/4OdQesAdyzV9oAXdRbYcqwDglqRsRFUuRKX3HZky906IGQVlBztdM11pqeuMP590hWcM1gxbrpiRQnSohT0FNcaVtRss8+5U28/+VxVgWfeAej3/bjCZnetjytVqkrW9Gf57DVQd94ytpzBa0F3k62NqUuey1rXqwJxb+36gA4LgwofV/qcPQ0OP0mQH1iiXTfp8NYLRaLoQbDFz3RyVjvf5L3O9a8zU5RA7RmVr/EMaVB2DxEkw6xbX+rn4/yD1dFVs41/nqugZnQTMMLSgu4CUkq+PVZIuSkgp3wTmIJhxU/8XnXaRihRoroFPftu1s85CFlOv85zRmmHNjXMzMAl4b08RZXUt3jPEHACX/0NFwwCExMKyJ1xPIhcYpiZgx5wHjRXw5nfg96nw1Plq/Ub+NuNtP4UYdoIupWRXXjXNbdYhv/ehkjqKapq5PfgTBBImXwlhcf1fJIRKUmQKgO3Pw9HP1PGjn6oq7CGxMOUaj9uuGZ6kxYayaEISbVbJS54IYXSFjPlw5ya47O+qcPXIaYPrJyQaVryuyiuOnKbCJ/O/Viusn1oEL69QAyCNyww7Qf/+yzu5/B9f8PGBkiG/9wd7iwmnkatMn6oDp/cyGdobCePgzB8BEl79FnzxKLxjj9s94x4I8mD5Mc2w5+YzMgH49+bjXhnIdCM6DWbepETZHUwmmLkSvvM53J+rBH7+PSoy7OAaeP4yaG00xORTiWEn6LMzVO6UV7YNfcKfD/eVcHvAWkKtdaoCjCu1GRf+DE67WCUsWvf/oOYEJM9UE6caTT+cMSaOScmRlNe38Np2P0x0FRKt8sdc+Du483MVOVO0E9b80NuWDTuGnaBfPj2ZQLOJjUfKKKxu6r9xfRlsfATW3qfCA53NrXJkHbz3Y3j/fhU3297KNyeqEMW7uT3AHpWy6AHXDDeZ4dr/wAW/g9EL1Uq9lW/ozIqaARFCcNdCtaz+8Q3Z3h+le5LY0XDDapU+Y/fLnVFg/XH0M/U3vvY+5crsj+oT6hvyh/+jMk66mkTPxxEeiW91gtmzZ8tt2wY3AXLPi9+wZncRdy0cw/1Lxvfe6NAHaqlxS23nsbGLVa6Vvlwc1jZYcy/seKH78ZBYdpLFqMY9RIlGtdBi2ZODsl2jGQxWm2Tpoxs5WFzHzy4az52eKIDhS3z1BHzwMwgfAXd/pQpz9ERKFT755aPdj0+4DC76I0R2qc/a1qzabXxEhU06yFigBloDzYX5EEKI7VLK2b2eG46CvuNEFcse/5LQQDMbf3ruyRWDDqyBV25SCa/GLMKWPAvb1qcIaK6kKmEOB89/nqjwCKJDLUSFWDCbBFarFcs7dxG4/zVkQDC2M36A2WxRCylK93V03TpqEYErXlIhiRrNEPLZ4TK+9czXBFtMrPnemYxNjDipzaHiOr7ILudoeT1NrTaiQixMSo7k3PGJxIY5H5EipURKMJmcqFTkCWw2ePYiyPtKxa5f8fjJbdY9CF/8VQUcnPE9td38uErlGxQJ5/5CLXoq2A4f/xoq7eX0JlymkpBtfVql902aAt96p//EeFLC8S/Ut/fKoyo3TVSqSp8weqFr82DFeyHhNLWCdhC4LehCiCXA3wAz8JSU8g89zgv7+YuBRuBmKWW/pUrcEXSAW5/byicHS1k2I4W/XNdl2X3hDnjmImhvonHuD3hUXs9/vz5BbEs+qwN/ywhRxXvW0/le2/exdXicJL8JeI6bAtZRL4NZ2fpz9ohxZMaHYbPZCKw4yBhRyJWLzmDReUucK8el0XiAH72ykze+KSA9NpR/3TSb00ZEUF7fwru7Cnl1Wz77i2p7vS4owMT1p6fznXNGMzIqpM/+d+VV89iGbD4/XIZNSianRHHNrDSumpVCUICTC4iMovwIPLEArC2w4jXIOr/z3OZ/qBXYpgBY/mLn4r6afHjvJyqHe08SxsPFf4JR9hw2dcXw3FKoyIYRU+CmPkS9cCe8/1OVBrg3zEEqf9PsW9WIvz99KPgGnr9UrT257gWVgsFF3BJ0IYQZOAycD+QDW4HrpZT7u7S5GPgeStDnAn+TUvaTrcd9Qc8pq+eSRzfR1Gbll0sncNtZo1US/qcWQV0RO+KWcn3pjTS3qfeXERfK6aFF/Lr8x4TKRt63XMBvuI3KJivfM73KPaY3aJEW7hE/Z5tpClWNnUv1I4MDuG/JeFbOyxi0vRqNEdS3tLN81Wb2FijhjgsLpKJLit2oEAuLJyQxKTmS8OAAyutb2JxTwcYjalFbYICJWxZk8t2FY4kK6Rwhbsut5NFPlJD3RmpMCD++YByXT0sZ2lH7pr/Cxw+qjI+3fQwxGaqYzBu3qfPLVsG0Hus4pIRDa5XbpvQARIxUkTmzbj45br62SIl6ZQ6MnK7qpzrcL01VsP639vUiEsISlLt15DQ1J1aeDUc+tMfO23U0fpwS9mnLT3YTnfgKXr5Bxd9Plll0cgAACIBJREFUvlplZTW5Po3prqDPB34lpbzQ/vrnAFLK33dp80/gUynlS/bXh4CFUsqivvodtKBXHlMLdJb+mVf31XHfa7sBWDw6mIer7iOxKYcttgnc2Ppz2ghg8YRE7j53bGdloWMb4YWr1H/95BkQGA65G1V+6Gv/o/7TAo2t7Rwta8Bqk4xNDCcsSNfT1vgGDS3tPPTeAd7ckU9zm42gABNnjInj6llpLJ6Y2OtI+kBRLX//5Ahr9xQDEB4UwAUTk0iMDGbLsYqOlNChgWZWzs/g2wtGERJo5pODpTz2STZHSusBGD8ighvmpjM1NZqgABNWm6SsvoWy2hZK65qpaGhVrhohEAICTAKzSRBgNnXs90VXLXLsCtnOZXu/T3r11zRY4iiJmMioyk0IJJ9n3suO1Bvd/jzDWkq5es8dRDfnUxs0gl0jr8Nia2Jq0auEtlVhFWZ2jVzOlvTbaQ042bUS0VzMpJK3mFTyFuGt6h9nuymIw/Hnkx81C5O0klb9NVnlH2PCRm7MGWyd9xjXzB3cPIi7gn41sERKeZv99UpgrpTyni5t1gB/kFJusr9eD9wvpexTsQct6M9fBsc+U+lpL/8Hb3yTz2/e2sHj8mHOMO8n25bMNW2/4owpWdy9cCwTkyNP7iP3C1Xd3FHQNjACLntULRTSaIYJzW1WapvaiA4NJDDAuZHerrxq/vD+QTYfreh2PCIogJsXZHLrglHE9PC1W22SN77J5y/rDlNYM/TL9CNo5JnAPzLHdBiAdmniz+1X87jVuPxHSVTyz8C/MN2U0+34Ftt4ftl2K0dk6oB9BNDOItM3rDCv52zznpPOt0kzT1sv4k/t1zI5LZ637l4wKFvdFfRrgAt7CPrpUsrvdWnzHvD7HoL+Uynl9h593QHcAZCenj7r+PFBJOcpOwRPnqVG2Of/FtLn0/7hLwnI/4oGSxyfnf0Ss6dPIzFiAN9Uc21n5ris83XqWs0pxZGSOjYfraCmsY1RCWGcNz6R0MD+v4U2t1n5YG8x6w6UcKysAZtdOxIigkiICCIxIpj48EBMQmCzT6q22yRWm82+lbTbuutNz/F6V/ez6HJWYCOr8lPCWivIj5xJRZjxUT4mWzvjKj4mrWYb7aYgjsacyfHoeYOaM4tuymN8+YdEN+cBUB46lsNxi6gNVjozIiqEG+amD8pO/3K5QGdIU1ciRsINr8DIqYPrU6PRaIYB/Qm6M47hrUCWEGIUUAAsB27o0eYd4B4hxMuoSdGa/sTcbebdBeFJ8MXflNsk6wI4536Vg1yj0WhOUQYUdClluxDiHuBDVNjiM1LKfUKIO+3nnwTWoiJcslFhiy7m1BwEk6/UPm+NRqPpglOhG1LKtSjR7nrsyS77ErjbWNM0Go1G4wrDLpeLRqPRaHpHC7pGo9H4CVrQNRqNxk/Qgq7RaDR+ghZ0jUaj8RO0oGs0Go2foAVdo9Fo/ASvFbgQQpQBg0jm4jHigXJvG+Ek2lbPoG31DNpWY8mQUib0dsJrgu5rCCG29ZUfwdfQtnoGbatn0LYOHdrlotFoNH6CFnSNRqPxE7Sgd7LK2wa4gLbVM2hbPYO2dYjQPnSNRqPxE/QIXaPRaPwELegajUbjJ5wSgi6EWCKEOCSEyBZC/KyX8/cJIXbaf/YKIaxCiFj7uVwhxB77uUHWzHPazmeEEKVCiL19nBdCiEft72O3EGJml3P9vkcv2LrCbuNuIcSXQohpXc4N2WfqpK0LhRA1XZ6BB7qc87XP1SeeVfv90oQQG4QQB4QQ+4QQP+iljdefWSft9Jnn1S2klH79g6qylAOMBgKBXcDEftpfCnzS5XUuED9Etp4NzAT29nH+YuB9VG3decCWwbzHIbL1DCDGvn+Rw9ah/kydtHUhsMbdZ2cobO3R1mvPqv1+I4GZ9v0I4HDPz8cXnlkn7fSZ59Wdn1NhhH46kC2lPCqlbAVeBi7vp/31wEtDYlkPpJSfA5X9NLkc+LdUfAVECyFG4vp79LitUsovpZRV9pdfAametKc/nPhc+8LnPtceeO1ZBZBSFkkpv7Hv1wEHgJQezbz+zDpjpy89r+5wKgh6CpDX5XU+Jz90AAghQoElwOtdDkvgIyHEdiHEHR6z0jn6ei9Ov0cv8W3UKM2BL32mDuYLIXYJId4XQkyyH/PZz9XXnlUhRCYwA9jS45RPPbP92NmV4fC89opTNUWHOaKXY33Fal4KfCGl7DpCWiClLBRCJALrhBAH7aMob9DXe3HlPQ4pQohzUX8gZ3Y57EufKcA3qPwY9UKIi4G3gCx8+HPFh55VIUQ46h/LvVLK2p6ne7nEK8/sAHY62gyH57VPToURej6Q1uV1KlDYR9vl9PgKK6UstG9LgTdRXxW9RV/vxZX3OGQIIaYCTwGXSykrHMd97DNFSlkrpay3768FLEKIeHz0c7XjE8+qEMKCEsn/Sinf6KWJTzyzTtg5bJ7X/jgVBH0rkCWEGCWECET9IbzTs5EQIgo4B3i7y7EwIUSEYx+4AOg1+mCIeAe4yR45MA+okVIW4eR7HEqEEOnAG8BKKeXhLsd97TNFCDFCCCHs+6ej/i4q8MHPFXznWbV/Zk8DB6SUf+6jmdefWWfsHE7Pa3/4vctFStkuhLgH+BA1s/6MlHKfEOJO+/kn7U2XAR9JKRu6XJ4EvGn/Ww8AXpRSfuApW4UQL6EiLuKFEPnAg4Cli51rUVED2UAjcEt/79FTdjpp6wNAHPC4/fNrlyqL3ZB+pk7aejVwlxCiHWgClksV3uCLnyv4wLNqZwGwEtgjhNhpP/YLIL2Lvb7wzDpjp888r+6gl/5rNBqNn3AquFw0Go3mlEALukaj0fgJWtA1Go3GT9CCrtFoNH6CFnSNRqPxE7SgazQajZ+gBV2j0Wj8hP8PdaWqlDN8NakAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hddm.analyze.plot_posterior_nodes(m.nodes_db.loc[[\"v(easy)\", \"v(hard)\"], \"node\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [--------------------------150%---------------------------] 3 of 2 complete in 5.9 sec"
     ]
    }
   ],
   "source": [
    "ppc_data = hddm.utils.post_pred_gen(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>observed</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>SEM</th>\n",
       "      <th>MSE</th>\n",
       "      <th>credible</th>\n",
       "      <th>quantile</th>\n",
       "      <th>mahalanobis</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>stat</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>accuracy</th>\n",
       "      <td>0.910000</td>\n",
       "      <td>0.927500</td>\n",
       "      <td>0.050637</td>\n",
       "      <td>0.000306</td>\n",
       "      <td>0.002870</td>\n",
       "      <td>True</td>\n",
       "      <td>31.500000</td>\n",
       "      <td>0.345593</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean_ub</th>\n",
       "      <td>0.935440</td>\n",
       "      <td>0.951484</td>\n",
       "      <td>0.098785</td>\n",
       "      <td>0.000257</td>\n",
       "      <td>0.010016</td>\n",
       "      <td>True</td>\n",
       "      <td>46.400002</td>\n",
       "      <td>0.162416</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std_ub</th>\n",
       "      <td>0.421473</td>\n",
       "      <td>0.464852</td>\n",
       "      <td>0.111991</td>\n",
       "      <td>0.001882</td>\n",
       "      <td>0.014424</td>\n",
       "      <td>True</td>\n",
       "      <td>37.700001</td>\n",
       "      <td>0.387345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10q_ub</th>\n",
       "      <td>0.501000</td>\n",
       "      <td>0.521309</td>\n",
       "      <td>0.039991</td>\n",
       "      <td>0.000412</td>\n",
       "      <td>0.002012</td>\n",
       "      <td>True</td>\n",
       "      <td>32.200001</td>\n",
       "      <td>0.507850</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30q_ub</th>\n",
       "      <td>0.686000</td>\n",
       "      <td>0.658764</td>\n",
       "      <td>0.057687</td>\n",
       "      <td>0.000742</td>\n",
       "      <td>0.004070</td>\n",
       "      <td>True</td>\n",
       "      <td>71.500000</td>\n",
       "      <td>0.472135</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50q_ub</th>\n",
       "      <td>0.832000</td>\n",
       "      <td>0.820630</td>\n",
       "      <td>0.088531</td>\n",
       "      <td>0.000129</td>\n",
       "      <td>0.007967</td>\n",
       "      <td>True</td>\n",
       "      <td>59.900002</td>\n",
       "      <td>0.128429</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70q_ub</th>\n",
       "      <td>1.008000</td>\n",
       "      <td>1.053636</td>\n",
       "      <td>0.134238</td>\n",
       "      <td>0.002083</td>\n",
       "      <td>0.020102</td>\n",
       "      <td>True</td>\n",
       "      <td>40.200001</td>\n",
       "      <td>0.339962</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>90q_ub</th>\n",
       "      <td>1.573000</td>\n",
       "      <td>1.542863</td>\n",
       "      <td>0.247801</td>\n",
       "      <td>0.000908</td>\n",
       "      <td>0.062314</td>\n",
       "      <td>True</td>\n",
       "      <td>59.599998</td>\n",
       "      <td>0.121618</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean_lb</th>\n",
       "      <td>-1.049667</td>\n",
       "      <td>-0.990737</td>\n",
       "      <td>0.350941</td>\n",
       "      <td>0.003473</td>\n",
       "      <td>0.126632</td>\n",
       "      <td>True</td>\n",
       "      <td>33.798283</td>\n",
       "      <td>0.167918</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std_lb</th>\n",
       "      <td>0.430255</td>\n",
       "      <td>0.297619</td>\n",
       "      <td>0.251784</td>\n",
       "      <td>0.017592</td>\n",
       "      <td>0.080988</td>\n",
       "      <td>True</td>\n",
       "      <td>74.248924</td>\n",
       "      <td>0.526784</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10q_lb</th>\n",
       "      <td>0.491400</td>\n",
       "      <td>0.726807</td>\n",
       "      <td>0.332695</td>\n",
       "      <td>0.055416</td>\n",
       "      <td>0.166103</td>\n",
       "      <td>True</td>\n",
       "      <td>9.442060</td>\n",
       "      <td>0.707575</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30q_lb</th>\n",
       "      <td>0.799600</td>\n",
       "      <td>0.819943</td>\n",
       "      <td>0.335312</td>\n",
       "      <td>0.000414</td>\n",
       "      <td>0.112848</td>\n",
       "      <td>True</td>\n",
       "      <td>61.266094</td>\n",
       "      <td>0.060669</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50q_lb</th>\n",
       "      <td>1.130000</td>\n",
       "      <td>0.928716</td>\n",
       "      <td>0.356847</td>\n",
       "      <td>0.040515</td>\n",
       "      <td>0.167855</td>\n",
       "      <td>True</td>\n",
       "      <td>80.686699</td>\n",
       "      <td>0.564064</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70q_lb</th>\n",
       "      <td>1.192800</td>\n",
       "      <td>1.076653</td>\n",
       "      <td>0.398996</td>\n",
       "      <td>0.013490</td>\n",
       "      <td>0.172688</td>\n",
       "      <td>True</td>\n",
       "      <td>69.849785</td>\n",
       "      <td>0.291099</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>90q_lb</th>\n",
       "      <td>1.516800</td>\n",
       "      <td>1.302383</td>\n",
       "      <td>0.520954</td>\n",
       "      <td>0.045975</td>\n",
       "      <td>0.317368</td>\n",
       "      <td>True</td>\n",
       "      <td>71.995705</td>\n",
       "      <td>0.411585</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          observed      mean       std       SEM       MSE credible  \\\n",
       "stat                                                                  \n",
       "accuracy  0.910000  0.927500  0.050637  0.000306  0.002870     True   \n",
       "mean_ub   0.935440  0.951484  0.098785  0.000257  0.010016     True   \n",
       "std_ub    0.421473  0.464852  0.111991  0.001882  0.014424     True   \n",
       "10q_ub    0.501000  0.521309  0.039991  0.000412  0.002012     True   \n",
       "30q_ub    0.686000  0.658764  0.057687  0.000742  0.004070     True   \n",
       "50q_ub    0.832000  0.820630  0.088531  0.000129  0.007967     True   \n",
       "70q_ub    1.008000  1.053636  0.134238  0.002083  0.020102     True   \n",
       "90q_ub    1.573000  1.542863  0.247801  0.000908  0.062314     True   \n",
       "mean_lb  -1.049667 -0.990737  0.350941  0.003473  0.126632     True   \n",
       "std_lb    0.430255  0.297619  0.251784  0.017592  0.080988     True   \n",
       "10q_lb    0.491400  0.726807  0.332695  0.055416  0.166103     True   \n",
       "30q_lb    0.799600  0.819943  0.335312  0.000414  0.112848     True   \n",
       "50q_lb    1.130000  0.928716  0.356847  0.040515  0.167855     True   \n",
       "70q_lb    1.192800  1.076653  0.398996  0.013490  0.172688     True   \n",
       "90q_lb    1.516800  1.302383  0.520954  0.045975  0.317368     True   \n",
       "\n",
       "           quantile  mahalanobis  \n",
       "stat                              \n",
       "accuracy  31.500000     0.345593  \n",
       "mean_ub   46.400002     0.162416  \n",
       "std_ub    37.700001     0.387345  \n",
       "10q_ub    32.200001     0.507850  \n",
       "30q_ub    71.500000     0.472135  \n",
       "50q_ub    59.900002     0.128429  \n",
       "70q_ub    40.200001     0.339962  \n",
       "90q_ub    59.599998     0.121618  \n",
       "mean_lb   33.798283     0.167918  \n",
       "std_lb    74.248924     0.526784  \n",
       "10q_lb     9.442060     0.707575  \n",
       "30q_lb    61.266094     0.060669  \n",
       "50q_lb    80.686699     0.564064  \n",
       "70q_lb    69.849785     0.291099  \n",
       "90q_lb    71.995705     0.411585  "
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hddm.utils.post_pred_stats(data, ppc_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The returned data structure is a pandas `DataFrame` object with a hierarchical index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>rt</th>\n",
       "      <th>response</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>node</th>\n",
       "      <th>sample</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"10\" valign=\"top\">wfpt(easy)</th>\n",
       "      <th rowspan=\"10\" valign=\"top\">0</th>\n",
       "      <th>0</th>\n",
       "      <td>0.481109</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.755106</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.713106</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.100101</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.905104</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.716106</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>-0.873104</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.404109</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.566114</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.419107</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           rt  response\n",
       "node       sample                      \n",
       "wfpt(easy) 0      0  0.481109         1\n",
       "                  1  0.755106         1\n",
       "                  2  0.713106         1\n",
       "                  3  1.100101         1\n",
       "                  4  0.905104         1\n",
       "                  5  0.716106         1\n",
       "                  6 -0.873104         0\n",
       "                  7  0.404109         1\n",
       "                  8  1.566114         1\n",
       "                  9  1.419107         1"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ppc_data.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first level of the `DataFrame` contains each observed node. In this case the easy condition. If we had multiple subjects we would get one for each subject.\n",
    "\n",
    "The second level contains the simulated data sets. Since we simulated 500, these will go from 0 to 499 -- each with generated from a different parameter value sampled from the posterior.\n",
    "\n",
    "The third level is the same index as used in the data and numbers each trial in your data.\n",
    "\n",
    "For more information on how to work with hierarchical indices, see the [Pandas documentation](http://pandas.pydata.org/pandas-docs/stable/indexing.html#hierarchical-indexing-multiindex).\n",
    "\n",
    "There are also some helpful options like `append_data` you can pass to `post_pred_gen()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function post_pred_gen in module kabuki.analyze:\n",
      "\n",
      "post_pred_gen(model, groupby=None, samples=500, append_data=False, add_model_parameters=False, progress_bar=True)\n",
      "    Run posterior predictive check on a model.\n",
      "    \n",
      "    :Arguments:\n",
      "        model : kabuki.Hierarchical\n",
      "            Kabuki model over which to compute the ppc on.\n",
      "    \n",
      "    :Optional:\n",
      "        samples : int\n",
      "            How many samples to generate for each node.\n",
      "        groupby : list\n",
      "            Alternative grouping of the data. If not supplied, uses splitting\n",
      "            of the model (as provided by depends_on).\n",
      "        append_data : bool (default=False)\n",
      "            Whether to append the observed data of each node to the replicatons.\n",
      "        progress_bar : bool (default=True)\n",
      "            Display progress bar\n",
      "    \n",
      "    :Returns:\n",
      "        Hierarchical pandas.DataFrame with multiple sampled RT data sets.\n",
      "        1st level: wfpt node\n",
      "        2nd level: posterior predictive sample\n",
      "        3rd level: original data index\n",
      "    \n",
      "    :See also:\n",
      "        post_pred_stats\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(hddm.utils.post_pred_gen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to compute the summary statistics over each simulated data set and compare that to the summary statistic of our actual data by calling `post_pred_stats()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ppc_compare = hddm.utils.post_pred_stats(data, ppc_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          observed      mean       std       SEM       MSE credible  \\\n",
      "stat                                                                  \n",
      "accuracy  0.910000  0.927500  0.050637  0.000306  0.002870     True   \n",
      "mean_ub   0.935440  0.951484  0.098785  0.000257  0.010016     True   \n",
      "std_ub    0.421473  0.464852  0.111991  0.001882  0.014424     True   \n",
      "10q_ub    0.501000  0.521309  0.039991  0.000412  0.002012     True   \n",
      "30q_ub    0.686000  0.658764  0.057687  0.000742  0.004070     True   \n",
      "50q_ub    0.832000  0.820630  0.088531  0.000129  0.007967     True   \n",
      "70q_ub    1.008000  1.053636  0.134238  0.002083  0.020102     True   \n",
      "90q_ub    1.573000  1.542863  0.247801  0.000908  0.062314     True   \n",
      "mean_lb  -1.049667 -0.990737  0.350941  0.003473  0.126632     True   \n",
      "std_lb    0.430255  0.297619  0.251784  0.017592  0.080988     True   \n",
      "10q_lb    0.491400  0.726807  0.332695  0.055416  0.166103     True   \n",
      "30q_lb    0.799600  0.819943  0.335312  0.000414  0.112848     True   \n",
      "50q_lb    1.130000  0.928716  0.356847  0.040515  0.167855     True   \n",
      "70q_lb    1.192800  1.076653  0.398996  0.013490  0.172688     True   \n",
      "90q_lb    1.516800  1.302383  0.520954  0.045975  0.317368     True   \n",
      "\n",
      "           quantile  mahalanobis  \n",
      "stat                              \n",
      "accuracy  31.500000     0.345593  \n",
      "mean_ub   46.400002     0.162416  \n",
      "std_ub    37.700001     0.387345  \n",
      "10q_ub    32.200001     0.507850  \n",
      "30q_ub    71.500000     0.472135  \n",
      "50q_ub    59.900002     0.128429  \n",
      "70q_ub    40.200001     0.339962  \n",
      "90q_ub    59.599998     0.121618  \n",
      "mean_lb   33.798283     0.167918  \n",
      "std_lb    74.248924     0.526784  \n",
      "10q_lb     9.442060     0.707575  \n",
      "30q_lb    61.266094     0.060669  \n",
      "50q_lb    80.686699     0.564064  \n",
      "70q_lb    69.849785     0.291099  \n",
      "90q_lb    71.995705     0.411585  \n"
     ]
    }
   ],
   "source": [
    "print(ppc_compare)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, we did not have to define the summary statistics as by default, `HDDM` already calculates a bunch of useful statistics for RT analysis such as the accuracy, mean RT of the upper and lower boundary (ub and lb respectively), standard deviation and quantiles. These are listed in the rows of the DataFrame.\n",
    "\n",
    "For each distribution of summary statistics there are multiple ways to compare them to the summary statistic obtained on the observerd data. These are listed in the columns. `observed` is just the value of the summary statistic of your data. `mean` is the mean of the summary statistics of the simulated data sets (they should be a good match if the model reproduces them). `std` is a measure of how much variation is produced in the summary statistic.\n",
    "\n",
    "The rest of the columns are measures of how far the summary statistic of the data is away from the summary statistics of the simulated data. `SEM` = standard error from the mean, `MSE` = mean-squared error, `credible` = in the 95% credible interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can also tell `post_pred_stats()` to return the summary statistics themselves by setting `call_compare=False`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ppc_stats = hddm.utils.post_pred_stats(data, ppc_data, call_compare=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                   accuracy   mean_ub    std_ub    10q_ub    30q_ub    50q_ub  \\\n",
      "node       sample                                                               \n",
      "wfpt(easy) 0           0.86  1.062849  0.724248  0.484509  0.667707  0.905104   \n",
      "           1           0.98  0.977835  0.474633  0.498661  0.680059  0.882057   \n",
      "           2           0.82  0.996223  0.515579  0.555150  0.671149  0.817147   \n",
      "           3           0.88  0.875679  0.305067  0.510363  0.733061  0.835860   \n",
      "           4           0.92  0.815038  0.503699  0.510257  0.593757  0.626257   \n",
      "\n",
      "                     70q_ub    90q_ub   mean_lb    std_lb    10q_lb    30q_lb  \\\n",
      "node       sample                                                               \n",
      "wfpt(easy) 0       1.078501  1.749323 -1.398971  0.610027  0.740706  1.084301   \n",
      "           1       1.133453  1.642467 -0.624060  0.000000  0.624060  0.624060   \n",
      "           2       1.004145  1.785162 -1.002371  0.407241  0.611750  0.756748   \n",
      "           3       0.951658  1.376258 -0.928192  0.198288  0.706861  0.757861   \n",
      "           4       0.779755  1.143250 -1.067014  0.635796  0.646857  0.666057   \n",
      "\n",
      "                     50q_lb    70q_lb    90q_lb  \n",
      "node       sample                                \n",
      "wfpt(easy) 0       1.268101  1.595516  2.214336  \n",
      "           1       0.624060  0.624060  0.624060  \n",
      "           2       0.802147  1.145143  1.375747  \n",
      "           3       0.916359  1.091857  1.161356  \n",
      "           4       0.733756  0.934758  1.753776  \n"
     ]
    }
   ],
   "source": [
    "print(ppc_stats.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `DataFrame` has a row for each simulated data set. The columns are the different summary statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using PPC for model comparison with the `groupby` argument\n",
    "\n",
    "One useful application of PPC is to perform model\n",
    "comparison. Specifically, you might estimate two models, one for which\n",
    "a certain parameter is split for a condition (say drift-rate ``v`` for\n",
    "hard and easy conditions to stay with our example above) and one in\n",
    "which those conditions are pooled and you only estimate one\n",
    "drift-rate.\n",
    "\n",
    "You then want to test which model explains the data better to assess\n",
    "whether the two conditions are really different. To do this, we can\n",
    "generate data from both models and see if the pooled model\n",
    "systematically misses aspects of the RT data of the two\n",
    "conditions. This is what the ``groupby`` keyword argument is\n",
    "for. Without it, if you ran ``post_pred_gen()`` on the pooled model\n",
    "you would get simulated RT data which was not split by\n",
    "conditions. Note that while the RT data will be split by condition,\n",
    "the exact same parameters are used to simulate data of the two\n",
    "conditions as the pooled model does not separate them. It simply\n",
    "allows us to match the two conditions present in the data to the\n",
    "jointly simulated data more easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_pooled = hddm.HDDM(data)  # v does not depend on conditions\n",
    "m_pooled.sample(1000, burn=20)\n",
    "ppc_data_pooled = hddm.utils.post_pred_gen(m_pooled, groupby=[\"condition\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You could then compare ``ppc_data_pooled`` to ``ppc_data`` above (by\n",
    "passing them to ``post_pred_stats``) and find that the model with\n",
    "separate drift-rates accounts for accuracy (``mean_ub``) in both\n",
    "conditions, while the pooled model can't account for accuracy in\n",
    "either condition (e.g. lower ``MSE``)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining your own summary statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also define your own summary statistics and pass them to `post_pred_stats()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ppc_stats = hddm.utils.post_pred_stats(\n",
    "    data, ppc_data, stats=lambda x: np.mean(x), call_compare=False\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>stat</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>node</th>\n",
       "      <th>sample</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">wfpt(easy)</th>\n",
       "      <th>0</th>\n",
       "      <td>0.718194</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.945797</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.636476</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.659214</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.664474</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       stat\n",
       "node       sample          \n",
       "wfpt(easy) 0       0.718194\n",
       "           1       0.945797\n",
       "           2       0.636476\n",
       "           3       0.659214\n",
       "           4       0.664474"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ppc_stats.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `stats` can also be a dictionary mapping the name of the summary statistic to its function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary statistics relating to outside variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another useful way to apply posterior predictive checks is if you have trial-by-trial measure (e.g. EEG brain measure). In that case the `append_data` keyword argument is useful.\n",
    "\n",
    "Lets add a dummy column to our data. This is going to be uncorrelated to anything but you'll get the idea."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy.random import randn\n",
    "\n",
    "data[\"trlbytrl\"] = randn(len(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No model attribute --> setting up standard HDDM\n",
      "Set model to ddm\n",
      " [-----------------100%-----------------] 1 of 1 complete in 0.0 sec1.4 sec"
     ]
    }
   ],
   "source": [
    "m_reg = hddm.HDDMRegressor(data, \"v ~ trlbytrl\")\n",
    "m_reg.sample(1000, burn=20)\n",
    "\n",
    "ppc_data = hddm.utils.post_pred_gen(m_reg, append_data=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats import linregress\n",
    "\n",
    "ppc_regression = []\n",
    "for (node, sample), sim_data in ppc_data.groupby(level=(0, 1)):\n",
    "    ppc_regression.append(\n",
    "        linregress(sim_data.trlbytrl, sim_data.rt_sampled)[0]\n",
    "    )  # slope\n",
    "\n",
    "orig_regression = linregress(data.trlbytrl, data.rt)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                rt_sampled  response_sampled  index     rt  response  \\\n",
      "node sample                                                            \n",
      "wfpt 0      0     1.121020                 1      0  0.934       1.0   \n",
      "            1     0.487028                 1      1  0.802       1.0   \n",
      "            2     1.383025                 1      2  1.394       1.0   \n",
      "            3     0.762025                 1      3  1.213       1.0   \n",
      "            4     1.609036                 1      4 -1.434       0.0   \n",
      "...                    ...               ...    ...    ...       ...   \n",
      "            95    0.717025                 1     95  1.015       1.0   \n",
      "            96    0.711025                 1     96  0.827       1.0   \n",
      "            97    0.805024                 1     97  0.468       1.0   \n",
      "            98    1.552033                 1     98 -0.612       0.0   \n",
      "            99    0.672026                 1     99  1.138       1.0   \n",
      "\n",
      "                subj_idx condition  trlbytrl  \n",
      "node sample                                   \n",
      "wfpt 0      0          0      easy  0.588377  \n",
      "            1          0      easy -0.247001  \n",
      "            2          0      easy -0.347119  \n",
      "            3          0      easy  2.098002  \n",
      "            4          0      easy -0.850838  \n",
      "...                  ...       ...       ...  \n",
      "            95         0      hard  2.381048  \n",
      "            96         0      hard  0.181995  \n",
      "            97         0      hard  0.374229  \n",
      "            98         0      hard  0.278482  \n",
      "            99         0      hard  1.971242  \n",
      "\n",
      "[100 rows x 8 columns]\n",
      "                rt_sampled  response_sampled  index     rt  response  \\\n",
      "node sample                                                            \n",
      "wfpt 1      0     0.859229                 1      0  0.934       1.0   \n",
      "            1     0.430233                 1      1  0.802       1.0   \n",
      "            2     0.728231                 1      2  1.394       1.0   \n",
      "            3     0.570233                 1      3  1.213       1.0   \n",
      "            4     0.584233                 1      4 -1.434       0.0   \n",
      "...                    ...               ...    ...    ...       ...   \n",
      "            95    2.195267                 1     95  1.015       1.0   \n",
      "            96   -0.518234                 0     96  0.827       1.0   \n",
      "            97    0.573233                 1     97  0.468       1.0   \n",
      "            98    0.705231                 1     98 -0.612       0.0   \n",
      "            99    0.475234                 1     99  1.138       1.0   \n",
      "\n",
      "                subj_idx condition  trlbytrl  \n",
      "node sample                                   \n",
      "wfpt 1      0          0      easy  0.588377  \n",
      "            1          0      easy -0.247001  \n",
      "            2          0      easy -0.347119  \n",
      "            3          0      easy  2.098002  \n",
      "            4          0      easy -0.850838  \n",
      "...                  ...       ...       ...  \n",
      "            95         0      hard  2.381048  \n",
      "            96         0      hard  0.181995  \n",
      "            97         0      hard  0.374229  \n",
      "            98         0      hard  0.278482  \n",
      "            99         0      hard  1.971242  \n",
      "\n",
      "[100 rows x 8 columns]\n",
      "                rt_sampled  response_sampled  index     rt  response  \\\n",
      "node sample                                                            \n",
      "wfpt 2      0     0.646321                 1      0  0.934       1.0   \n",
      "            1     0.688321                 1      1  0.802       1.0   \n",
      "            2     0.417322                 1      2  1.394       1.0   \n",
      "            3    -0.748320                 0      3  1.213       1.0   \n",
      "            4     0.575322                 1      4 -1.434       0.0   \n",
      "...                    ...               ...    ...    ...       ...   \n",
      "            95    0.631321                 1     95  1.015       1.0   \n",
      "            96    0.563322                 1     96  0.827       1.0   \n",
      "            97    0.737320                 1     97  0.468       1.0   \n",
      "            98    0.640321                 1     98 -0.612       0.0   \n",
      "            99    0.515322                 1     99  1.138       1.0   \n",
      "\n",
      "                subj_idx condition  trlbytrl  \n",
      "node sample                                   \n",
      "wfpt 2      0          0      easy  0.588377  \n",
      "            1          0      easy -0.247001  \n",
      "            2          0      easy -0.347119  \n",
      "            3          0      easy  2.098002  \n",
      "            4          0      easy -0.850838  \n",
      "...                  ...       ...       ...  \n",
      "            95         0      hard  2.381048  \n",
      "            96         0      hard  0.181995  \n",
      "            97         0      hard  0.374229  \n",
      "            98         0      hard  0.278482  \n",
      "            99         0      hard  1.971242  \n",
      "\n",
      "[100 rows x 8 columns]\n"
     ]
    }
   ],
   "source": [
    "cnt = 0\n",
    "for (node, sample), sim_data in ppc_data.groupby(level=(0, 1)):\n",
    "    print(sim_data)\n",
    "    cnt += 1\n",
    "    if cnt > 2:\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'slope')"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAQaklEQVR4nO3de6xlZXnH8e+vM3gDjUM54Milg82kFqkXOFqq1tqCKYJ1+KOkULUTSzKxsV6amnaoTUn/IB2qNbVJNZ0AdRoJlCoNk+JtOpVaY0FnALk4yoBQmDJljhdQNEHBp3/sNfF4OOPZZ1/OnvOe7yeZ7LXetdZez5Mz8zvvrL332qkqJElt+ZlJFyBJGj3DXZIaZLhLUoMMd0lqkOEuSQ1aPekCAI455phat27dpMuQtFi7d/94+fTTJ1fHCrV79+5vVNXUfNsOi3Bft24du3btmnQZkhYr+fGy/4aXXJL/OdQ2L8tIUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDotPqEp6qnWbb5jYue/fcu7Ezq3RcOYuSQ0y3CWpQYa7JDXIcJekBhnuktSgBd8tk+RK4A3Agao6tRt7H/BbwA+Ae4G3VtUj3baLgYuAJ4F3VtWnx1S7tCQm+a4VaVD9zNw/Apw9Z2wHcGpVvRi4G7gYIMkpwAXAi7pjPpRk1ciqlST1ZcFwr6rPAd+aM/aZqnqiW70JOKFb3gBcU1WPV9V9wD3AK0ZYrySpD6O45v77wCe75eOBB2dt29eNSZKW0FDhnuS9wBPAVQeH5tmtDnHspiS7kuyamZkZpgxJ0hwDh3uSjfReaH1TVR0M8H3AibN2OwF4aL7jq2prVU1X1fTU1Lxf3i1JGtBA4Z7kbOBPgTdW1fdnbdoOXJDk6UlOBtYDXxy+TEnSYvTzVsirgdcCxyTZB1xC790xTwd2JAG4qareVlV3JbkW+Aq9yzVvr6onx1W8JGl+C4Z7VV04z/AVP2X/S4FLhylKkjQcP6EqSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq0ILhnuTKJAeS3Dlr7OgkO5Ls7R7XzNp2cZJ7knwtyW+Oq3BJ0qH1M3P/CHD2nLHNwM6qWg/s7NZJcgpwAfCi7pgPJVk1smolSX1ZMNyr6nPAt+YMbwC2dcvbgPNmjV9TVY9X1X3APcArRlSrJKlPg15zP66q9gN0j8d248cDD87ab1839hRJNiXZlWTXzMzMgGVIkuYz6hdUM89YzbdjVW2tqumqmp6amhpxGZK0sg0a7g8nWQvQPR7oxvcBJ87a7wTgocHLkyQNYtBw3w5s7JY3AtfPGr8gydOTnAysB744XImSpMVavdAOSa4GXgsck2QfcAmwBbg2yUXAA8D5AFV1V5Jrga8ATwBvr6onx1S7JOkQFgz3qrrwEJvOPMT+lwKXDlOUJGk4fkJVkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoKHCPckfJbkryZ1Jrk7yjCRHJ9mRZG/3uGZUxUqS+jNwuCc5HngnMF1VpwKrgAuAzcDOqloP7OzWJUlLaNjLMquBZyZZDTwLeAjYAGzrtm8DzhvyHJKkRRo43Kvqf4H3Aw8A+4FHq+ozwHFVtb/bZz9w7HzHJ9mUZFeSXTMzM4OWIUmaxzCXZdbQm6WfDDwfODLJm/s9vqq2VtV0VU1PTU0NWoYkaR7DXJY5C7ivqmaq6ofAdcArgYeTrAXoHg8MX6YkaTGGCfcHgDOSPCtJgDOBPcB2YGO3z0bg+uFKlCQt1upBD6yqm5N8DLgFeAK4FdgKHAVcm+Qier8Azh9FoZKWzrrNN/S13/0DHPNTn2/LuUM/h3oGDneAqroEuGTO8OP0ZvGSpAnxE6qS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDhrq3jLRURnFTKmklceYuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYNFe5JnpvkY0m+mmRPkl9JcnSSHUn2do9rRlWsJKk/w87cPwh8qqpeCLwE2ANsBnZW1XpgZ7cuSVpCA4d7kucArwGuAKiqH1TVI8AGYFu32zbgvGGLlCQtzjAz9xcAM8A/Jrk1yeVJjgSOq6r9AN3jsfMdnGRTkl1Jds3MzAxRhiRprmHCfTVwGvDhqnoZ8D0WcQmmqrZW1XRVTU9NTQ1RhiRprmHCfR+wr6pu7tY/Ri/sH06yFqB7PDBciZKkxRo43Kvq/4AHk/xCN3Qm8BVgO7CxG9sIXD9UhZKkRRv2O1TfAVyV5GnA14G30vuFcW2Si4AHgPOHPIckaZGGCvequg2YnmfTmcM8ryRpOH5CVZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KChwz3JqiS3Jvm3bv3oJDuS7O0e1wxfpiRpMUYxc38XsGfW+mZgZ1WtB3Z265KkJTRUuCc5ATgXuHzW8AZgW7e8DThvmHNIkhZv9ZDH/y3wJ8CzZ40dV1X7Aapqf5Jj5zswySZgE8BJJ500ZBlaKus23zDpEiT1YeCZe5I3AAeqavcgx1fV1qqarqrpqampQcuQJM1jmJn7q4A3JjkHeAbwnCQfBR5Osrabta8FDoyiUElS/waeuVfVxVV1QlWtAy4A/qOq3gxsBzZ2u20Erh+6SknSoozjfe5bgNcl2Qu8rluXJC2hYV9QBaCqbgRu7Ja/CZw5iueVtLJM6gX7+7ecO5HzjpOfUJWkBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVo4HBPcmKSzybZk+SuJO/qxo9OsiPJ3u5xzejKlST1Y5iZ+xPAH1fVLwJnAG9PcgqwGdhZVeuBnd26JGkJDRzuVbW/qm7plr8L7AGOBzYA27rdtgHnDVukJGlxRnLNPck64GXAzcBxVbUfer8AgGMPccymJLuS7JqZmRlFGZKkztDhnuQo4OPAu6vqO/0eV1Vbq2q6qqanpqaGLUOSNMtQ4Z7kCHrBflVVXdcNP5xkbbd9LXBguBIlSYs1zLtlAlwB7KmqD8zatB3Y2C1vBK4fvDxJ0iBWD3Hsq4C3AHckua0b+zNgC3BtkouAB4DzhytRkrRYA4d7VX0eyCE2nzno80qShucnVCWpQYa7JDXIcJekBhnuktQgw12SGjTMWyE1Ies23zDpEiQd5py5S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAb5VkhJK94k3158/5Zzx/K8ztwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBfkJ1CH5phqTDlTN3SWrQ2GbuSc4GPgisAi6vqi3jOpczaEn6SWOZuSdZBfw98HrgFODCJKeM41ySpKca12WZVwD3VNXXq+oHwDXAhjGdS5I0x7guyxwPPDhrfR/wy7N3SLIJ2NStPpbka4s8xzHANwaucHlaaT2vtH5hmfWc2SuXvWGQp1hW/Y7IT/Scy4Z6rp871IZxhXvmGaufWKnaCmwd+ATJrqqaHvT45Wil9bzS+oWV1/NK6xeWrudxXZbZB5w4a/0E4KExnUuSNMe4wv1LwPokJyd5GnABsH1M55IkzTGWyzJV9USSPwQ+Te+tkFdW1V0jPs3Al3SWsZXW80rrF1ZezyutX1iinlNVC+8lSVpW/ISqJDXIcJekBi2bcE9ydJIdSfZ2j2vm2ecZSb6Y5MtJ7kryl5OodVT67PnEJJ9Nsqfr+V2TqHUU+um32+/KJAeS3LnUNY5CkrOTfC3JPUk2z7M9Sf6u2357ktMmUeco9dHzC5P8d5LHk7xnEjWOUh/9vqn72d6e5AtJXjLyIqpqWfwB/hrY3C1vBi6bZ58AR3XLRwA3A2dMuvYx97wWOK1bfjZwN3DKpGsfV7/dttcApwF3TrrmAXpcBdwLvAB4GvDluT8v4Bzgk93f5zOAmydd9xL0fCzwcuBS4D2TrnkJ+n0lsKZbfv04fsbLZuZO7/YF27rlbcB5c3eonse61SO6P8v5FeN+et5fVbd0y98F9tD7hPBytGC/AFX1OeBbS1XUiPVza44NwD91f59vAp6bZO1SFzpCC/ZcVQeq6kvADydR4Ij10+8Xqurb3epN9D4LNFLLKdyPq6r90As0er/pnyLJqiS3AQeAHVV18xLWOGp99XxQknXAy+j9j2U5WlS/y9R8t+aY+8u4n32Wk9b6Wchi+72I3v/URuqw+rKOJP8OPG+eTe/t9zmq6kngpUmeC/xrklOr6rC9NjuKnrvnOQr4OPDuqvrOKGobh1H1u4wteGuOPvdZTlrrZyF995vk1+mF+6tHXcRhFe5VddahtiV5OMnaqtrf/Rf1wALP9UiSG4GzgcM23EfRc5Ij6AX7VVV13ZhKHYlR/oyXqX5uzdHa7Tta62chffWb5MXA5cDrq+qboy5iOV2W2Q5s7JY3AtfP3SHJVDdjJ8kzgbOAry5ZhaPXT88BrgD2VNUHlrC2cViw3wb0c2uO7cDvde+aOQN49ODlqmVqpd2OZMF+k5wEXAe8paruHksVk35leRGvQP8ssBPY2z0e3Y0/H/hEt/xi4Fbgdnqz9b+YdN1L0POr6f2X73bgtu7POZOufVz9dutXA/vpvfi2D7ho0rUvss9z6L2r6V7gvd3Y24C3dcuh92U39wJ3ANOTrnkJen5e97P8DvBIt/ycSdc9xn4vB74969/srlHX4O0HJKlBy+myjCSpT4a7JDXIcJekBhnuktQgw12SGmS4a8VLcmOSFfUlzWqf4S5JDTLctaIkOTLJDd09/+9M8jtztl+Y5I5u22Wzxh9L8jdJbkmyM8lUN/7zST6VZHeS/0rywqXuSZqP4a6V5mzgoap6SVWdCnzq4IYkzwcuA34DeCnw8iQHbzt8JHBLVZ0G/CdwSTe+FXhHVZ0OvAf40NK0If10hrtWmjuAs5JcluRXq+rRWdteDtxYVTNV9QRwFb0vBgH4EfDP3fJHgVd3d+J8JfAv3W2m/4Hel6dIE3dY3RVSGrequjvJ6fTu/fFXST4za/N8t2o95FPRmxw9UlUvHWWN0ig4c9eK0l16+X5VfRR4P72v6zvoZuDXkhyTZBVwIb1LMND7t/Lb3fLvAp+v3n3z70tyfvfcGct3YUoDcOauleaXgPcl+RG9u0r+Ab2Qp3r3kb8Y+Cy9WfwnqurgbYe/B7woyW7gUeDgC7FvAj6c5M/pfa3jNfS+M1OaKO8KKfUhyWNVddSk65D65WUZSWqQM3dJapAzd0lqkOEuSQ0y3CWpQYa7JDXIcJekBv0/ZGyBzoYP6QQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(ppc_regression)\n",
    "plt.axvline(orig_regression, c=\"r\", lw=3)\n",
    "plt.xlabel(\"slope\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the simulated data sets have on average no correlation to our trial-by-trial measure (just as in the data) but we also get a nice sense of the uncertainty in our estimation."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "ff3096d2709bbb36a4584c44f6e6ffdb5e175071e94d34047f50b078bfdc1c6d"
  },
  "kernelspec": {
   "display_name": "hddmnn_tutorial",
   "language": "python",
   "name": "hddmnn_tutorial"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
